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SENSITIVITY OF ITERATIVE SPECTRALLY SMOOTH 
TEMPERATURE/EMISSIVITY SEPARATION TO INSTRUMENT NOISE 



5 RELATED APPLICATION 

This application claims the benefit of U.S. provisional application Serial No. 
60/221,238, filed July 27, 2000, entitled Sensitivity of Iterative Spectrally Smooth 
Temperature/Emissivity Separation to Instrument Noise. 
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TECHNICAL FIELD OF THE INVENTION 

This invention relates to determination of the error statistics for retrieved 
surface temperature and emissivity, and more particularly to determination of error 
statistics for retrieved surface temperature and emissivity by second-order analytical 
error propagation and Monte Carlo simulation. 
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BACKGROUND OF THE INVENTION 

The Iterative Spectrally Smooth Temperature-Emissivity Separation algorithm 
(ISSTES) is an inversion method for retrieving surface temperature and emissivity 
from remotely sensed hyperspectral thermal infrared radiance. It is important to 
quantify the error sensitivity for ISSTES - this is needed by system designers in 
determining system parameters required for a given retrieval error, and is also needed 
by consumers of the temperature/emissivity products as an indicator of product 
quality. 
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SUMMARY OF THE INVENTION 

A second-order analytical error model and Monte Carlo simulation are 
developed in detail and the results presented and compared. In one embodiment, the 
error statistics are computed for the case of a low-emissivity (0.65) surface at a 
5 temperature of 300 K, imaged by the SEBASS LWIR instrument from 2500 ft. 
through a warm humid atmosphere. With a noise level of 0.6 uW/cm 2 -sr-um 
(equivalent to an SNR of about 1200:1 in this scenario), the contribution of instrument 
noise to the standard deviation of retrieved temperature was 0.18 K, and the 
contribution to the standard deviation of retrieved spectral emissivity varied in the 

10 range 0.002 - 0.004 as a function of band number. The performance of ISSTES for 
sensors with lower SNR is discussed. ISSTES should give approximately unbiased 
surface temperature and emissivity estimates for noise levels below about 1 uW/cm 2 - 
sr-|^m (SNR > 750:1). Higher noise levels could result in biases due to non-linearity 
of the inversion within the radiance error bars. ISSTES converged successfully for 

15 noise levels up to about 5 uW/cm 2 -sr-um (SNR > 150:1), but showed significant 
probability of non-convergence at higher noise levels. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

For a more complete understanding of the present invention and the 
advantages thereof, reference is now made to the following description taken in 
connection with the accompanying drawings: 
5 FIGURE 1 is a plot of emissivity versus wavelength (microns) for a test panel 

emissivity; 

FIGURE 2 is a plot of radiance versus wavelength (microns) illustrating 
entrance aperture spectral radiance; 

FIGURE 3 is a flow chart of Monte Carlo simulation in accordance with the 
10 present invention; 

FIGURE 4 is a plot of mean temperature versus ensemble size (thousands) for 
Monte Carlo simulation; 

FIGURE 5 is a plot of temperature Sigma versus noise level illustrating 
standard deviation of estimated surface temperature; 
15 FIGURE 6 is a plot of temperature bias versus noise level illustrates surface 

temperature bias from analytical model evaluation showing the bias of the estimated 
surface temperature; 

FIGURE 7 is a three-dimensional plot of an analytical estimate of emissivity 
Sigma showing standard deviation of the estimate surface emissivity; 
20 FIGURE 8 is a three-dimensional plot for an analytical estimate of emissivity 

bias showing bias of the estimated surface emissivity; 

FIGURES 9a, 9b and 9c are plots of frequency versus temperature illustrating 
distribution of retrieved surface temperature for three scenarios; 

FIGURE 10 is a plot of temperature Sigma versus noise level for Monte Carlo 
25 simulation showing standard deviation of surface emissivity; 

FIGURE 1 1 is plot of temperature bias versus noise level illustrating surface 
temperature bias for Monte Carlo simulation; 

FIGURE 12 is a three-dimensional plot illustrating standard deviation of 
surface emissivity from Monte Carlo simulation; 
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FIGURE 13 is a three-dimensional plot showing bias of emissivity from 
Monte Carlo simulation; 

FIGURE 14 is a plot of temperature Sigma versus noise level illustrating the 
difference in analytical/Monte Carlo estimates of standard deviation; 
5 FIGURE 15 is a plot of temperature bias versus noise level illustrating the 

difference in analytical/Monte Carlo estimates of bias; 

FIGURE 16 is a three-dimensional plot showing differences in 
analytical/Monte Carlo estimates of emissivity standard deviation; 

FIGURE 17 is a three-dimensional plot illustrating differences in 
1 0 analytical/Monte Carlo estimates of emissivity bias; and 

FIGURE 18 is a plot of probability versus noise level illustrating frequency of 
convergence versus noise level. 
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DETAILED DESCRIPTION OF THE INVENTION 

An important problem in remote sensing is the retrieval of land surface 

temperature and spectral emissivity from airborne/spaceborne hyperspectral thermal 

infrared (TIR) data. A terrestrial surface with temperature t and emissivity s radiates 

5 a TIR spectral signature that propagates through the atmosphere to the instrument's 

entrance aperture, where it appears as the radiance L. (A symbol in boldface type 

denotes a vector or matrix quantity.) The forward model at the individual band level is 

given by: 

'V 

L= ' (1.1a) 
A, 

10 

L t = TfrB^i) + r,.(l - Si )LT + Zf (Lib) 
In equation Lib, i is band number (?' = l,..,N), B t (i) is the average Planck 
blackbody radiance in band i , and % i , L" p and Ej° wn are the atmospheric transmission 
15 between ground and sensor, atmospheric upwelling radiance at the sensor, and the 
atmospheric downwelling radiance plus solar radiance at the ground. 

The ISSTES algorithm provides an inversion method that maps L back to 
surface temperature and emissivity, denoted (t,e). (A circumflex over a symbol 
denotes a retrieved value.) Specifically, by iteration the value t that minimizes the 
20 spectral smoothness measure S : 
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(1.2a) 



L, -L? -r t Lf 



(1.2b) 



The retrieved emissivity e is then given by s = s(f ) . 



For most surface materials the error in the retrieved temperature and 



5 emissivity is small if L is known exactly. (It is also assumed that t , L" p and L down 
are known exactly). However, in practice L is not known exactly - there is at least 
some measurement noise, whose magnitude depends on the sensor design and 
operation. These errors in L propagate through the retrieval process and emerge as 
errors in estimated surface temperature and emissivity. 

10 In accordance with this invention, there is developed a method for estimating 

the errors in retrieved temperature and emissivity from the magnitude and distribution 
of the measurement noise. Such method is useful to instrument designers as a guide in 
determining the appropriate instrument noise level required for a given retrieval 
accuracy, and is also useful to consumers of the temperature/emissivity products as an 

1 5 indicator of product quality. 

The error statistics for retrieved surface temperature and emissivity is 
estimated in two ways: second-order analytical error propagation and Monte Carlo 
simulation. Each technique has advantages and disadvantages. The analytical model 
provides limited information about the distribution of temperature and emissivity 

20 errors, namely the first two moments. The estimate of the mean includes second-order 
terms in the Taylor series expansion of the retrieval function and can therefore test for 
biased retrievals. The estimate of the standard deviation includes first-order terms. By 
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postulating a structure for the covariance of the instrument noise and assuming that 
the retrieval function is quadratic within these error bars, the analytical model has the 
advantage of real-time performance and therefore can be used in an operational 
retrieval algorithm to provide product quality estimates. Monte Carlo simulation 

5 involves retrieving surface temperature and emissivity from sufficiently large 
ensembles of radiance spectra. This technique has the drawback of non-real-time 
performance, but yields more information than the analytical model, namely an 
empirical distribution of the errors. This distribution includes information not 
accounted for by the analytical model, such as the presence of statistical outliers 

10 resulting from non-convergence or incorrect convergence of the algorithm. 
Comparison of the two methods reinforces any conclusions implied by one method 
and identifies parameter ranges where the linear and quadratic approximations, 
inherent in analytical error propagation, begin to break down. 

In accordance with one embodiment, the simulation consisted of a surface at 

15 an assumed temperature of 300 K and with the low spectral emissivity as shown in 
Figure 1. Most other temperature/emissivity separation algorithms would be expected 
to perform poorly on this test case due to their underlying assumptions, whereas 
ISSTES should perform well. 

The terms r ,V P and l/ own of equation 1.2b were generated using 

20 MODTRAN with atmospheric profiles measured in central Oklahoma during June 
1997, with surface air temperature at about 300 K and relative humidity consistently 
near 80%. The resulting entrance aperture radiance is shown in Figure 2. 
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The second order analytical error model for ISSTES and the Monte Carlo error 
simulation technique will be discussed in detail. Monte Carlo simulation and 
analytical propagation techniques are applied to evaluate the error in retrieved surface 
temperature and emissivity arising from random instrument noise. 
5 Analytical Error Model for ISSTES, given a variable X (an iVxl vector) and a 

scalar-valued non-linear function y = F(x) , analytical error propagation provides a 
technique for estimating error statistics for y in terms of those for X and the partial 
derivatives of ^(x) . When the range of errors in X are small enough that F(x) can 
be accurately represented by a first-order approximation, F(x) is said to be "linear 
10 within the error bars for X", and it is only necessary to expand the function 
y = F(x) in a first order Taylor series about X : 

y + Sy*F(x)+— SX 2.1a 
dF 

Sy^—SX 2.1b 
y dX 

2 dF dF 1 

(7i« — — 2.ic 

* dX ffl dX 

15 

dF 

where — is the Jacobian matrix (lxiV) of the first partial derivatives 
dX 

ofF(x), <5X represents the perturbing random variable, represents the covariance 
matrix (NxN) of SX, <5y represents the resulting random error in the dependent 
variable, cr\ is the variance of Sy and superscript-t denotes matrix transpose. (The 
20 convention of expressing the derivative of a scalar function, F(x), with respect to a 
column vector, X , as a row vector has been adopted so that the Taylor series may be 
expressed without the use of the transpose operator.) Note that equation 2.1c is valid 
without restriction on the error distribution of SX , and that an unbiased noise source 
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(i.e., (<SX) = 0) results in an unbiased random error in the dependent variable, 

W = o. 

When the errors SX are so large that f(x) is not adequately represented by a 

first order expansion, F(x) is said to be "non-linear within the error bars for X ." In 

5 such cases a second order expansion can be undertaken: 

dF 1 , d 2 F 
dX 2 dX 2 

0* =&)*\ t™ce[ri ^ j 2.2b 

d 2 F 

where is the Hessian matrix ( N x N ) of the second partial derivatives of 

dX 

F(x) . Equation 2.2b assumes that errors in X have a Gaussian distribution and that 
1 0 (<5X) = 0 . The result is clear for uncorrected errors, where the matrix is diagonal. 

The result for correlated errors follows from the uncorrelated case by changing 

variables to a coordinate system in which T^is diagonal, and noting that the trace 

operation is invariant under the orthogonal transformation used in the diagonalization. 

Note that in general, a function F(x) that is non-linear within the error bars for X 
1 5 results in a biased random error in the dependent variable, (Sy) ^ 0 , even for an 

unbiased noise source (i.e., (SX) = 0). This resulting bias, denoted jS^, is given by 

the expression on the right hand side of equation 2.2b. 

If the function is not given explicitly in the form y = ^(x) , but instead is 

defined implicitly as the level set of another function G : 
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g(x, v )=o (2-3) 

Then under certain conditions on G (the implicit function theorem), equation 

2.3 implies that>> can be expressed as a function of X . However, it is usually difficult 

to write down an explicit representation of this function in closed form y = F(x) 

amenable to differentiation. Nevertheless, equations 2.1 and 2.2 can still be applied 

because what these equations require is not F(x) but its derivatives, and these can be 

written in terms of those for G(X, y). The first derivative of F(x)can be obtained by 

implicit differentiation of equation 2.3: 

dG dG dF A „ A ^ 

— + = 0 => (2.4a) 

dX dy dX 



dF^__(dG_\ dG 
dX ~ { dy ) 8X 



(2.4b) 



Note in equation 2.4b that is a scalar, and and are 1 x N vectors. 
dy dX dX 

Similarly, the second derivative of F(x) can be written by differentiating 

15 equation 2.4b: 



d 2 F = fdG\ 8 2 G | d 2 G 
dX 2 {dy)\dX 2+ 8ydX 



1 dF^f dGY 2 fdGVf d 2 G ^ 8 2 G dF} 
dx\ + [dy) {dXj [dXdy + dy 2 dX) (>) 



d 2 G d' 1 G d 2 G 

Note in equation 2.5 that — — is a scalar, and are 1 x N vectors, 

dy 2 dXdy dydX 

, d 2 F J d 2 G Ar „ . 
and and are N x N symmetric matrices. 

dx 2 dx 2 y 

This general approach will now be applied to estimate the covariance and bias 
of surface temperature and emissivity errors retrieved using ISSTES. Let the symbols 
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(i,e) denote both the inverse mapping implied by the algorithm and the retrieved 
values: 

M x 

L (t,e) (2.6) 

Recall that i is defined as the temperature that minimizes the smoothness 
measure S in equation 1.2. Near a solution, this condition can be re-expressed in 
terms of the derivative of S : 

ds(hj) n 



(2.7) 



dt 

Equation 2.7 is of the form in equation 2.3 and implicitly defines i as a 
function of L, i = f(L). Therefore equations 2.4 and 2.5 can be applied to ^ to write 

the derivatives — and : 
dL dL 2 

dt ( d 2 sY d 2 S 

^ = \W) aL¥ (2 - 8a) 



d L t = d z S 
dL 2 dt 2 



d 3 S ^ ( d 3 S V dt 



Y di ( d 2 s) z f d 2 s X( d 3 s d 3 s di) ,„ 

J [^J[i[^ + ^^:J (2 - 8b) 



In computing equations 2.8, it is helpful to write S in terms of the 
(N-2)x.N matrix Kof second order differences, so that (up to a multiplicative 
constant which can be ignored): 



(2.9a) 
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K = 



-1 2 -1 

-1 2 -1 



(2.9b) 



-1 2 -lj 

Differentiating equation 2.9 produces, via matrix operations, an analytically 
dS 



convenient expression for 



dt 



dS 



(2.10) 



where s and — are obtained from equation 1 .2. The higher-order derivatives 
dt 

of required in equations 2.8a, 2.8b can be readily written down in similar manner. 

The derivatives in equations 2.8a, 2.8b can then be substituted in equations 2.1 and 
1 0 2.2 to estimate the covariance and bias of i . 

The covariance and bias of retrieved emissivity can be estimated in a similar 

manner using equations 2.1 and 2.2 once and ^ £ ' have been derived. But 

dL dL 2 dL 

can be obtained by differentiating the definition e = s (t ) with respect to L : 



ds i _ ds t de; di 
dL dL dt dL 



(2.11) 



Note in equation 2. 1 1 that is a scalar, and , and — are 1 x N 
dt dL dL dL 



d 2 i. 

Finally, can be obtained by differentiating equation 2.1 1 : 
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d 2 s. t _d 2 £i 

dL 2 ~ dL 2 



ra 2 g,.V di (diX d 2 £i 

[dtdL I dL\dL j dLdt 



* * (2 , 2) 

5? <*L 2 dt 2 [dL) dL 



Note in equation 2.12 that is zero from equation 1.2, ^ £ ' and — are 

dL 2 dL 

lxN vectors, and ^ £ J and are NxN symmetric matrices. 
dL dL 2 y 

In summary, equations 2. 1 and 2.2 can be used to estimate the covariance and 

5 bias for surface temperature and emissivity retrieved by the ISSTES algorithm, where 

the required derivatives of i are given by equation 2.8, and the derivatives of e are 

given by equations 2.1 1 and 2.12. 

Given a variable X and a function y = F(x) as described previously, Monte 

Carlo simulation provides a technique for generating an empirical approximation of 

10 the error distribution for y in terms of the distribution for X and the function ^(x). 

The Monte Carlo simulation process is illustrated in Figure 3. 

The covariance of the measurement variables is used to generate an 

unbiased ensemble of measurement errors. These errors are added to the vector X , to 

produce an ensemble simulating the results of making the same measurement many 

15 times, each time with slightly different errors. This ensemble of measurements is then 

evaluated using the function F(x) to give an ensemble of dependent variables, which 

can then be used to estimate the mean, bias and variance of the dependent variable 

according to the following equations: 

5^=^|> (3.1a) 
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1 N 



(3.1b) 



h = y N -y 



(3.1c) 



This general approach will be applied to compute the error distribution for 
simulated measurements of a low emissivity material as described. Specifically, X is 
the radiance spectrum of Figure 2, and F(x) is the ISSTES retrieval function 
previously outlined, implemented using Newton's method. Each band (i.e. element of 
X) has an independent associated noise variable, and the generation of each element 
in an ensemble will involve selecting different pseudo-random numbers as the noise 
values for each band. The surface temperature and emissivity are retrieved from each 
spectrum in an ensemble, the temperature and emissivity errors are calculated, the 
error distribution is compiled, and the error statistics estimated according to equations 



The ensemble size N needs to be chosen large enough that statistically 
significant estimates of the surface temperature bias can be made. The specific 
condition on N is shown in equation 3.2 where j3 is the bias of equation 3.1c, e.g. 
for retrieved surface temperature. 



This was done informally by computing y N for increasing values of TV , and 
then plotting the difference y N -y^, An example is shown in Figure 4 for two noise 
levels. The choice of JV = 10,000 appears to satisfy the condition in equation 3.2 for 
surface temperature biases as low as 0.03 K. 



3.1. 



\y N -y\«P 



(3.2) 
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The following describes the results of applying analytical error and Monte 
Carlo simulation to the test case first discussed, measured in the presence of random 
instrument noise. It is assumed that the measured radiance L differs from the true 
radiance L by additive Gaussian noise, which is independent in each band. Denoting 
the noise in channel i by the random variable v t : 

L i =L i +v i i = l,....N (4.1) 
For simplicity it is also assumed that the noise level has the same value v in 
all bands, so that each v i has zero mean and variance v 2 . 

In the analytical error model, the condition above translates into a diagonal 
covariance matrix for the radiance error that has the following form: 

(4.2) 

The analytical estimates for the standard deviation and bias of surface 
1 5 temperature are shown in Figures 5 and 6. 

For a spectrometer such as SEBASS [2], the noise level v = 0.6 uW/cm 2 -sr- 
um gives an SNR of about 1200:1 in this example. This noise level results in a 
standard deviation of surface temperature of about 0.18 K and a bias of about 0.03 K. 
For other sensors, noise levels can be higher. At noise levels of 3 uW/cm 2 -sr- 
20 um, corresponding to an SNR of 240: 1, the standard deviation of surface temperature 
reaches about 0.88 K and the surface temperature bias reaches a value of about 0.70 
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K. The implication for such instruments, at least in this scenario, is that instrument 
noise levels should be kept below about 1 uW/cm 2 -sr-u.m if the ISSTES algorithm is 
to give approximately unbiased temperature retrievals. This corresponds to an SNR of 
720: 1 and implies a lower bound on the sensor integration time or spectral resolution. 
5 The analytical estimates for the standard deviation and bias of surface 

emissivity are shown in Figures 7 and 8. 

For SEBASS, the contribution to the standard deviation of surface emissivity 
is in the range 0.002 - 0.004, increasing near the edges of the wavelength region due 
to the loss of signal from the surface caused by increased atmospheric absorption at 
10 these wavelengths. The surface emissivity bias has the opposite sign as the 
temperature bias (as expected by conservation of energy), and is in the range -0.0004 
- -0.0008. 

For noise levels of 3 uW/cm 2 -sr-urn, standard deviation of surface emissivity 
reaches the 0.015-0.025 range, and the emissivity bias reaches values in the -0.009 - - 
15 0.015 range. As in the case of retrieved surface temperature, instrument noise levels 
should be kept below about 1 uW/cm 2 -sr-um if the ISSTES algorithm is to give 
approximately unbiased emissivity retrievals in this scenario. 

Radiance spectra ensemble sizes of 10,000 were generated for Monte Carlo 
simulation for different values of v, and the temperature and emissivity were retrieved 
20 using the ISSTES algorithm. Histograms for the retrieved surface temperature for v 
= 1,2 and 3 ^/cm 2 -sr-um are shown in Figures 9a, 9b and 9c. 

Figure 10 illustrates that the surface temperature standard deviation increases 
for increasing noise levels. Also, since the true surface temperature was 300 K, a bias 
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appears as the noise level increases. The surface temperature standard deviation and 
bias were estimated from these distributions according to equations 3.1 and in Figures 
10 and 11. 

From equations 2.1 and 2.2, standard deviation of surface temperature is 
5 approximately linear in noise level, and surface temperature bias is approximately 
quadratic in noise level. Good agreement with the analytical model was obtained even 
for relatively high noise levels. 

Standard deviation and bias of surface emissivity were computed from the 
ensembles as well, and plotted in Figures 12 and 13. The linear and quadratic 
10 dependence of the standard deviation and bias on noise levels is again evident from 
these figures, and excellent agreement with the analytical results was obtained. 

The differences between the analytical and Monte Carlo estimates for standard 
deviation and bias of retrieved surface temperature are shown in Figures 14 and 15, 
drawn at the same scale as the previous plots. These Figures indicate that the first- 
15 order analytical model for standard deviation and the second-order analytical model 
for bias are adequate for estimating error statistics of retrieved surface temperature, 
even for relatively high noise levels of 3 jaW/cm 2 -sr-urn. 

The differences between the analytical and Monte Carlo estimates for standard 
deviation and bias of retrieved surface emissivity are shown in Figures 16 and 17, 
20 drawn at the same scale as the previous plots. As expected from the case for surface 
temperature, the first-order analytical model for standard deviation and the second- 
order analytical model for bias are adequate for estimating error statistics of retrieved 
surface emissivity out to noise levels of 3 uW/cm 2 -sr-um. 
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As the noise level increases, the noise pattern begins to interfere with the 
spectral pattern of atmospheric features, so that at some noise value the algorithm will 
start failing to converge, or will converge but to an improbable value. The standard 
deviation and bias for retrieved temperature/emissivity for a given noise level were 
5 estimated from the analytical model, and a range of five standard deviations about the 
bias was chosen. If the retrieved surface temperature for a point in an ensemble was 
within this range, then the algorithm was deemed to have converged successfully. The 
frequency of successful convergence for ensembles with varying noise levels was 
calculated, and the results are shown in Figure 18. ISSTES showed successful 
10 convergence for noise levels up to about 5 |aW/cm 2 -sr-um, or equivalent SNR of 
about 150:1. 

A second-order analytical model and Monte Carlo simulation of the sensitivity 
of the ISSTES algorithm to uncorrelated gaussian instrument noise have been 
described. For a given scenario, these techniques estimate the standard deviation and 

1 5 bias of surface temperature and emissivity retrieved by the ISSTES algorithm, and the 
noise levels at which algorithm convergence begins to fail. As an example, error 
statistics were computed for the case of a low-emissivity (0.65) surface at 300 K 
imaged by the SEBASS LWTR instrument from 2500 ft. through a warm humid 
atmosphere. With an SNR of about 1200:1 in this case, the contribution of instrument 

20 noise to the standard deviation and bias of retrieved surface temperature were 0.18 K 
and 0.03 K respectively, and the contribution to the standard deviation and bias of 
retrieved surface emissivity were in the ranges 0.002 - 0.004 and -0.0004 - -0.0008 
respectively. 
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The standard deviation of retrieved surface temperature and emissivity scale is 
approximately linear up to a noise level of 3 |iW/cm 2 -sr-urn, or equivalently down to 
an SNR of about 240:1 for this scenario. The bias also varies approximately 
quadratically with noise level within the same range. ISSTES should give 

5 approximately unbiased surface temperature and emissivity estimates for noise levels 
below about 1 uW/cm 2 -sr-um (SNR > 720:1) Higher noise levels could result in 
significant biases, due to non-linearity of the inversion within the radiance error bars. 

Based on the close agreement between the analytical and Monte Carlo error 
estimates for the ISSTES algorithm, the analytical error computation can be integrated 

10 into an ISSTES software package. This software not only performs surface 
temperature and emissivity retrievals on SEBASS data but also supplies the standard 
deviation of retrieved values as a product quality figure-of-merit. 

While the invention has been described in connection with a preferred 
embodiment, it is not intended to limit the scope of the invention to the particular 

15 form set forth, but, on the contrary, it is intended to cover alternatives, modifications, 
equivalents as may be included within the spirit and scope of the invention as defined 
in the appended claims. 



